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1.  Introduction 


This  report  considers  the  problem  of  multi-target  localization  using  transient  signals  from  a 
single-sensor  as  well  as  multi-sensor  point  of  view.  Here  it  is  assumed  that  target  identification  is 
not  possible,  and  therefore,  no  association  between  measurements  and  targets  are  available. 
Furthermore,  the  number  of  targets  in  the  surveillance  region  is  unknown.  Additionally,  due  to 
the  limited  range  of  the  sensors,  missed  detections  can  occur  and  the  presence  of  clutter  can 
induce  false  alarms.  An  example  of  such  a  scenario  is  shooter  localization  using  a  network  of 
acoustic  gunfire  detection  systems  (GDSs)  (1 ).  The  individual  GDSs  composed  of  a  passive 
array  of  microphones  are  able  to  localize  a  gunfire  event  by  measuring  the  direction  of  arrival  for 
both  the  acoustic  wave  generated  by  the  muzzle  blast  and  the  shockwave  generated  by  the 
supersonic  bullet  (2-5).  Due  to  echo,  reverberation,  and  the  dissipative  nature  of  the  acoustic 
signal,  missed  detections  and  false  alarms  are  prevalent  in  acoustic  source  localization. 
Furthermore,  due  to  the  transient  nature  of  the  event,  continuous  observations  are  not  available 
and  recursive  Bayesian  tracking  schemes  cannot  be  employed. 

Conventional  multi-target  tracking  (MTT)  approaches  like  Multiple  Hypothesis  Tracking 
(MHT)  (6,  7)  or  the  Joint  Probability  Data  Association  (JPDA)  filter  (8,  9),  which  address  the 
data-association  problem,  either  are  too  computationally  demanding  or  cannot  be  applied  to  the 
transient  event  localization  problem  since  these  approaches  require  persistent  measurement 
signals  and  a  fixed  number  of  targets.  In  MHT,  all  possible  combinations  of  tracks  and  data 
associations  are  exhaustively  evaluated;  therefore,  it  is  an  impractical  scheme  since  the  number  of 
mappings  between  data  and  targets  will  grow  exponentially  with  the  number  of  targets  (JO,  11). 
Though  more  efficient  than  MHT,  JPDA  methods  are  not  optimal  since  the  detection  is  performed 
separately  from  tracking  and  cannot  initiate  tracks  at  low  signal-to-clutter  ratios  (12).  A 
multi-target  extension  of  the  simultaneous  localization  and  mapping  (SLAM)  problem  for 
streaming  data  is  presented  by  Garcia-Femandez  et  al.  (13).  The  multi-target  simultaneous 
localization  and  mapping  (MSLAM)  scheme  is  based  on  the  parallel  partition  particle  filter  and  it 
outperforms  the  well-known  FastSLAM  (14)  when  there  are  multiple  targets  in  the  surveillance 
area.  Jensfelt  and  Kristensen  (15)  discuss  an  extension  of  the  MHT,  known  as  the 
multi-hypothesis  localization  (MHL),  for  mobile  robot  localization.  The  MHL  uses  a 
multi-hypothesis  Kalman  filter  along  with  a  probabilistic  formulation  of  hypothesis  correctness  to 
generate  and  track  Gaussian  hypotheses.  However,  almost  all  of  the  MTT  techniques  are 
recursive  algorithms  that  require  persistent  observations  and  are  futile  in  dealing  with  transient 
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signals. 


A  finite  point  process  approach  known  as  the  probability  hypothesis  density  (PHD)  filter, 
introduced  by  Mahler  (16),  allows  a  more  tractable  implementation  of  multi-target  tracking 
approaches  since  it  only  propagates  the  first-order  moment  of  the  multi-target  density.  Moreover, 
the  PHD  filter  is  able  to  handle  a  time-varying  number  of  targets,  missed  detections,  and  false 
alarms.  Though  the  track  labeling  problem  is  not  considered  in  the  PHD  filter,  a  track  labeling 
method  combined  with  the  PHD  approach,  is  proposed  by  Lin  et  al.  (17)  for  multi-target  tracking. 
Since  the  implementation  of  an  exact  PHD  filter  is  intractable,  a  sequential  Monte  Carlo  (SMC) 
or  particle  filtering  approach  (18)  and  the  Gaussian  sum  filtering  scheme  (19)  have  been  devised 
to  approximate  the  PHD  filter.  Convergence  properties  for  the  particle  PHD  filter  and  Gaussian 
mixture  PHD  filter  are  presented  in  references  (20)  and  (21),  respectively.  An  SMC 
implementation  of  a  finite  set  statistical  filter  for  the  localization  of  an  unknown  number  of 
speakers  in  a  multipath  environment  using  time  difference  of  arrival  (TDOA)  measurements  has 
also  been  proposed  (18,  22-24).  Similar  to  traditional  MTT  algorithms,  the  PHD  filter  is  a 
recursive  Bayesian  approach,  which  also  requires  persistent  observations  for  track  update.  A 
finite  point  process  approach  to  maximum  likelihood  based  multi-target  localization  of  an 
unknown  number  of  targets  from  transient  signals  has  not  been  considered  yet. 

Traditionally,  multi-target  localization  involves  the  maximum  likelihood  based  approach,  where 
the  selected  model  yields  the  maximum  likelihood  of  observing  the  given  data  across  a  possible 
number  of  targets  and  all  possible  target-data  associations.  As  the  number  of  sensors  increases, 
the  possible  combination  of  target-data  associations  dramatically  increases,  and  the  problem  often 
becomes  intractable.  Development  of  a  multi-target  detection  and  localization  scheme  based  on  a 
probabilistic  framework  known  as  modeling  field  theory  (MFT)  is  presented  by  Deming  and 
Perlovsky  (25).  Though  the  computational  complexity  of  a  MFT-based  approach  scales  linearly 
with  the  size  of  the  problem,  it  involves  an  iterative  scheme  similar  to  the  expectation 
maximization  (EM)  and  an  ad-hoc  likelihood  ratio  test  is  needed  to  prune  the  number  of  targets. 
An  iterative  maximum  likelihood  optimization  technique  based  on  a  modified  deterministic 
annealing  EM  (MDAEM)  algorithm  for  multi-target  localization  and  velocity  estimation  using 
TDOAs  is  given  by  Carevic  (26).  Since  the  MDAEM  algorithm  is  executed  for  an  assumed 
number  of  targets,  Carevic  (27)  provides  a  systematic  approach  for  determining  which  of  the 
target  models  estimated  by  the  MDAEM  algorithm  are  related  to  the  true  targets.  Both  the 
MFT-based  approach  and  the  MDAEM  algorithm  require  that  the  assumed  number  of  targets  is 
greater  than  or  equal  to  the  true  number  of  targets.  Also,  the  measurements  are  only  assumed  to 
contain  clutter/false  alarms  and  the  problem  of  missed  detection  is  not  considered. 
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For  the  multi-target  localization  problem  considered  here,  we  use  the  frequentist  counterpart  to 
the  Bayesian  filtering  approach,  i.e.,  the  maximum  likelihood  algorithm.  The  localization 
problem  is  formulated  in  two  dimensions  and  the  measurements  considered  here  are  the  range  and 
the  bearing  to  the  targets.  Each  sensor  acquires  several  range  and  bearing  measurement  pairs  and 
the  proposed  algorithm  estimates  the  number  of  targets  and  their  corresponding  locations  based 
on  the  erroneous  measurements.  The  number  of  targets  is  certainly  different  from  the  number  of 
measurements  due  to  clutter  and  missed  detection.  The  proposed  approach  is  a  twofold  scheme 
that  includes  an  EM  algorithm  to  estimate  the  target  locations  for  a  given  number  of  targets  and 
an  information  theoretic  algorithm  to  select  the  target  model,  i.e.,  number  of  targets  and  their 
locations.  The  main  advantage  of  the  proposed  scheme  is  that  it  scales  linearly  with  the  size  of 
the  problem  and  avoids  the  curse  of  dimensionality  associated  with  the  traditional  MHT-based 
multi-target  localization  scheme.  For  example,  if  there  are  N  targets,  S  sensors,  and  m 
measurements  per  sensor,  the  computational  complexity  for  the  proposed  scheme  is  in  the  order 
of  O  ( NSm )  while  the  computational  complexity  for  the  traditional  scheme  is  in  the  order  of 
O  (iV5m) .  Unlike  the  methods  used  by  Deming  and  Perlovsky  (25)  and  Carevic  (26),  the 
proposed  approach  accounts  for  probability  of  detection  and  missed  detection  along  with  clutter 
and  false  alarms. 

The  structure  of  this  report  is  as  follows:  formulation  of  multi-target  localization  problem  and 
the  corresponding  solution  for  the  single-sensor  scenario  are  presented  in  sections  2,  and  3, 
respectively.  A  numerical  example  demonstrating  the  single-sensor  algorithm  is  presented  in 
section  4.  Formulation  of  multi-target  localization  problem  and  the  corresponding  solution  for 
the  multi-sensor  scenario  is  presented  in  sections  5,  and  6,  respectively.  Numerical  simulation 
demonstrating  the  multi-sensor  algorithm  are  presented  in  section  7.  Section  8  presents  the 
results  obtained  from  implementing  the  proposed  algorithm  on  experimental  data.  Finally, 
Section  9  concludes  the  report  and  discusses  the  current  research  challenges. 


2.  Problem  Formulation 


Consider  a  two-dimensional  (2-D)  scenario  where  there  are  N  targets  located  in  R2.  The  target 
locations  are  denoted  as 


(Ti,  T2,  . . . ,  T^) 


T 

±  21 

T 

±yi 


T 


22 


T, 


V2 


(1) 


Due  to  clutter  and  miss  detection,  the  number  of  targets  and  the  target  location  set  that  may 
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induce  a  measurement  are  jointly  modeled  as  a  finite  point  process  known  as  the  Poisson  Point 
Process  (PPP)  (28).  A  realization  of  this  PPP,  S,  is  denotes  as 


*  = 


^X2 

ty2 


txm 

tym 


(m,{ ti,  t2,  tm}). 


(2) 


The  points  in  the  PPP,  i.e.,  t,  occur  in  a  state  space  S  =  M2  and  the  realizations  of  the  PPP  is 
defined  in  a  bounded  subset  7Z  C  S.  It  is  assumed  that  T  =  (Ti,  T2,  . . . ,  TN)  E  7Z.  In  a  PPP, 
the  points  t.,  are  not  necessarily  distinct  and  their  order  is  irrelevant.  The  PPP  is  fully 
parameterized  by  the  intensity  function: 


N 

A(t)  =  J^pD(t)(5(t  -  T^, 

i= 1 


(3) 


where  S(-)  is  the  Dirac  delta  function,  and  pD(-)  represents  the  probability  of  detection.  The 
number  of  points  in  the  PPP,  S,  is  determined  by  sampling  the  discrete  Poisson  random  variable 
M  with  probability  mass  function  (pmf): 


Pu(m)  =  (A exP  {-  A(t)dt}  . 


(4) 


The  m  points,  tj  E  1Z,  j  —  1, . . . ,  m,  are  defined  as  i.i.d.  samples  of  a  random  variable  T  on  1Z 
with  probability  density  function  (pdf): 


Pt  (f ) 


m 

JnX(t)dt 


(5) 


The  number  of  targets  and  their  locations  are  unknown,  but  a  sensor  located  at 


SX  Sy 


is  able 


to  observe  the  range  and  bearing  to  the  targets.  The  measurement  equation  is,  Vj  =  1, . . . ,  m, 


2  j  =  h(tXj,tVj)  +  wj, 


'  y/it*,  -  sxy  +  (tw  -  Sy)2 

arctan ((tyj  -  Sy)/(tXj  -  Sx)) 


Wr. 

+ 

3 

(6) 


The  noise  w j  is  assumed  to  be  i.i.d.  samples  of  Gaussian,  zero  mean  process  with  variance 

k2  (f 


£w  = 


0  a2 


.  Thus,  the  likelihood  can  be  written  as 

/(z|t)  —J\f(  z  |  h(t),  Ew) . 


(7) 
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Let’s  assume  that  the  each  point  t,-  in  the  PPP  realization  £  is  observed  by  the  sensor  and 
subsequently,  the  sensor  generates  a  measurement  zj  G  T  C  M2.  Let  -0  =  (tn,  {z1;  ....  zm}), 
then  if)  is  a  realization  of  a  PPP,  \P, defined  on  the  space  T,  with  intensity 

JV 

z/(z)  =  /  Z(z|t)A(t)dt  =  y~'pD(Tj)A/’(z  I  h(Tj),  Ew)  +  i/d(z),  (8) 

*=i 

where  /vc;(z)  is  the  clutter  intensity.  Thus  the  measured  points  are  modeled  as  either  samples 
from  one  of  the  N  Gaussian  densities  or  from  the  clutter  intensity.  Now  the  localization  problem 
can  be  formally  defined  as  follows: 

Given  =  (m,  { z , ,  . . . ,  zm}),  a  realization  of  the  PPP,  VP,  and  the  clutter  intensity  zy,/(z),  find 
an  estimate  of  the  number  of  targets,  N,  target  locations,  Ti,  T2,  •  •  • ,  Tn,  and  the 
corresponding  probability  of  detection,  pr)(T,). 


3.  Proposed  Single-Sensor  Solution 


The  proposed  twofold  solution  to  the  problem  defined  in  the  previous  section  consist  of  an  EM 
algorithm  to  estimate  the  target  locations  for  a  given  model  for  the  intensity  and  a  penalized  log 
likelihood  based  information  criteria  to  select  the  appropriate  model.  Here,  the  only  difference  in 
the  intensity  model  is  the  number  of  Gaussian  components  in  each  model  and  their  corresponding 
locations.  Details  of  the  twofold  solution  is  given  next. 

3.1  EM  Method 

Assume  the  model  is  known,  i.e.,  the  number  of  targets  is  given.  Given  the  number  of  targets,  N, 
the  intensity  z/(z  |  T)  is  the  superposition  of  N  weighted  Gaussian  components  along  with  the 
clutter  intensity: 

N 

I  T)  =  5>d(T,)  Af(z  I  h(Tj),  £w)  +  ud{z  I  T),  (9) 

i— 1 

where  the  parameter  vector  of  the  i-th  component  is  T,  and  T  =  ('T ,  ,  T2,  . . . ,  Tn)- 

3.1.1  E  -  Step 

The  natural  choice  of  the  “missing  data”  are  the  conditionally  independent  random  indices  kj, 
kj  G  {0,1,2,...,  iV },  that  identify  which  of  the  Gaussian  components  generated  the 
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measurement  zy  Here  k:J  =  0  indicates  that  the  measurement  is  generated  by  clutter.  Let 


zc  =  {m,  (z !,  fci), . . . ,  (zm,  fcm)} 


(10) 


denote  the  complete  data.  For  i  e  (0, 1,  2, . . . ,  N},  let 


zc(*)  =  {(zj,kj)  :  kj  =  i}. 


(ID 


Let  nc(i)  >  0  denote  the  number  of  indices  j  such  that  kj  =  i,  and  let 


ipc(i)  =  (wc(*)>  Zc(0  } 


(12) 


It  follows  from  the  definition  of  kj  that  ipc(i)  is  a  realization  of  the  PPP  whose  intensity  is 


^cCO 


PD( Ti)  Af(z\  h(Tj),  Ew ) ,  if  i  G  (1,  2, ... ,  iV}; 
Vd{  z|T),  if  z  =  0. 


(13) 


Now  only  considering  the  first  scenario,  i  e  {1,  2, . . . ,  N},  the  pdf  of  ipc(i)  can  be  written  as 


p(^c(i)  |  T)  =  exp  J  pD  (Ti)  AT  (z  |  h(T*),  Ew)dz  J  ]^[  pD(Ti)  J\f  (  Zj  |  h(Tj),  Ew ) 

j:kj=i 

(14) 

and 

pC0c(0)  |  T)  =  exp  J  isci(z  |  T)dz  ^  z/d(zj  |  T)  (15) 


J':fcj=0 


Since  the  superposed  components  are  independent,  we  have 


P  (zc  |  T)  =  p  (^c(0)  |  T)  p  (^c(l)  |  T) p  (^(2)  |  T) . . .  p  tyc(N)  |  T) 

=  exp  (-  f  u(z\T)dz\  pD(Tfc.)A/'(zi  |  h(Tfc.),  Ew  )  z/d(zj  |  T) 


j:kj=D 


The  log  likelihood  function  of  T  given  zc  is 


£(T|zc)  =  -  f  v(z\T)dz  +  log  PD( Tfc.)  +  Y  log  Af  ( z,  |  h(Tfc.),  Sw ) 

j'-kj^ti)  j'-kj^Q  (16) 

+  Y  log  ^(zj  i T) 

j:kj=$ 
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Note  that  the  conditional  pdf  of  the  missing  data,  (ki, ...  ,km),  given  T  and  z  can  be  written  as 


P  (ki, . . .  ,km\z.T) 


gjzcj  T) 
P  (z  I  T) 


tt  PD(TfcJ  )  J\f  ( z j  |  h(TjCj),  Sw  )  -|-|- 

J-J-  1S(Z  ;  I  T)  J-J- 


|  T) 
z/(Zi|T) 


Invoking  the  Poisson  Gambit,  p  (ki, . . . ,  km  |  z,  T)  can  be  written  as 


p(k1,...,km\z,T)=p(k1\z,T)p(k2\z,T)...p  (. km  |  z.  T) 


(17) 


where  the  individual  pdfs  can  be  written  as 


P(kj  |  z.T) 


pD(Tt.)Ar(z3-|h(Tt.),Ew) 
v(*S  I  T) 

v, d  (zj  |  T) 
u(z,  I  T)  > 


if  kj  ^  0; 
else. 


(18) 


Let  n  —  0,1,...  denote  the  EM  iteration  index,  and  let  the  current  feasible  value  of  T  be 
T(")  —  . . . ,  V  Now  the  EM  auxiliary  function  is  the  conditional  expectation 


Q  (T  |  T(n)) 


E 


ki  fcj7 


z,T(n) 


[£  (T 


N  N 


E  •••  E£ITI^) 

k\  — 0  km—fy 


n 


^(Tfcj)AT(Zj|h(TA,.),  Sw) 
is{zj  |  T) 


n 


^d(zj  I  T) 
u[zj  |  T) 


(19) 


Substituting  the  log  likelihood  yields 

Q  (T  |  T<">)  =  -  f  „(z|T)<iz  +  £ 

Jr  j= 1  fcl=l 

N 

■■■  jy  log  N  (  Zj  I  h(Tfc.),  Ew  ) 


m  N 


km  —  1 


V  (  z,  |  h 

V  (z j  I  TW) 


(20) 
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Note  that 


N 


Y  log  Af(zj  |h(Tfe.),  Ew) 

k\  ? 5 kj — i  j kj _|_ i ,  •  •  • , km — 1 

log  Af(zj  I  h(Tfc.),  Ew) 

Thus,  Q  (T  |  T^n^)  can  be  written  as 

N 

Q  (T  |  T(n))  =Y,Qi  (T*  |  T(n)) 


Ar(z<|h(T<”>'),E, 


V  z 


TW) 


■V  I  Z,  I  h  (Tgj  ■  Sw 

"(Zj  I™) 


2=1 


where 


(21) 


(22) 


Qi  (Ti  I  TW)  =  -  /  AT  ( z  I  h(Tj),  Ew  )  dz 


'r 


E 

3= 1 


-  AA(Zi|h(T;n)l,Ew 


z/  (zj  |  Thd) 


log  J\f  ( Zj  |  h(T?;),  Ew  ) 


(23) 


This  completes  the  E-step. 

3.1.2  M-Step 

The  M-step  maximizes  Q  (T  |  T("))  over  all  feasible  T,  i.e., 

T(n+1)  =  argrnax  Q  (T  |  T(ri))  .  (24) 

Assuming  there  is  no  functional  relation  between  T,  and  T j  for  i  ^  j,  the  required  M-step 
maximum  is  found  by  maximizing  the  expressions  Qi  (T,;  |  T(") )  separately.  Let 

T-n+1)  =  argrnax  Qt  (T,:  |  T(n))  ,  1  <  i  <  N.  (25) 


Before  we  further  proceed,  note  J\f  ( Zj  |  h(Tz),  Ew  )  can  be  written  as 


A/"(zj  |  h(Tj),  Ew  )  = 


\J  det(27rE^ 


:  exp 


_2  (ZJ  ”  h(Ti))TE“1  (zj  —  h(Tj)) 


(26) 
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Thus 


log U  ( Zj  I  h(Tj),  £w )  =  -  log  \J det(27rSw)  -  A  (Zj  -  h(T,))T  S”1  (Zj  -  h(T,)) .  (27) 


2  vr? 


Define  the  weight  cu*  (z j  |  Tlnl,  £w)  as 


to,  (Zj|T^,£w) 


TV  (z,|h  T  n)  ,£ 


y  i i  j  ^ w 

^  (zj  |  T(n)) 


The  weight  c u;  (z j  |  T,r9,  £w)  is  the  probability  that  the  point  z;  is  generated  by  the  i-th 
component  given  the  current  estimates  T*"-).  Now  equation  25  can  be  rewritten  as 


rp(n+l) 

1 1  —  arg  max 


P  m 

-  /  TV(z|h(T,),  £w)  dz  +  J]a;i(zj|T^,£w)logTV(zi|h(Ti),  £v 

JT  .7  =  1 


Then  T)n  satisfies  the  necessary  condition: 


Y  ^  (zi  I  TW>  Sw)  Swx  [Zj  -  h(Tj)]  VTih  (T?:) 


TV(z|h(T,),  Sw)^-1  [z  —  h(Tj)]  VTih  (Tj)  dz 


The  above  condition  may  be  simplified  to 


m  „ 

Y^(zJ\T(-n\^)[zJ-h(Tl)}  =  /  AT(z|h(T4),  £w)[z-h(Ti)]  dz  (31) 

3  =  1  '  T 

Based  on  the  assumptions  fr  J\f  (  z  |  h(Tj),  £w  )  z  dz  «  h(Tj)  and 

Jr  TV  (  z  |  h(Tj),  £w  )  dz  «  1,  the  EM  updates  T-n+1')  is  calculated  as  the  solution  to  the  equation 


X^(z,-|TW  £w)  z 

h  (T”+1)  «  ^ - 

X^(Zj|TW,£w) 
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Define 


y*  = 


Vn 

y*i 


m 

E  Ui  (zi  I  T(n)>  sw)  Z j 

3= 1 


m 


E  (zi 

3= 1 


TW,E, 


Now  Tjn+1) 


is  the  solution  to  the  nonlinear  equation 


\J (Xj"+1*  -  S+  +  (7’'”+1)  -  Syf 

Vri 

_arctan((Tir+1)  -  SX)/(T<)?+1)  -  Sy))_ 

y<i>i_ 

Thus  we  have 


rji(n+ 1) 

1  Xi 


f  y}i 

Y  1+tan2  (y^.  f 


1/2 


2/r, 


1+tan2  (^ .  ) 


+  'S'a:) 
1/2 

+  'S'xj 


if  -  f  <  y<+  <  f ; 

else, 


and 


(33) 


(34) 


(35) 


?J"+1)  =  tan  (j,*)  (rj”+1>  -  Sj  +  5„.  (36) 

Given  the  likelihood  function  is  bounded  above,  the  EM  iteration  nearly  always  converges  to  a 
local  maximum  of  the  likelihood  function.  However,  the  convergence  rate  is  only  linear  and  the 
first  few  iterations  are  commonly  observed  to  result  in  significant  improvements  in  the  likelihood 
function. 

3.2  Information  Criteria  for  Model  Selection 

The  previous  subsection  provides  an  approach  to  estimate  the  target  locations  given  the  number  of 
targets.  Presented  in  this  section  is  an  information  theoretic  approach  to  select  the  appropriate 
model,  i.e.,  the  number  of  targets  or  the  number  of  Gaussian  components.  Model  selection  can 
be  approached  in  terms  of  the  Kullback-Leibler  information  of  the  true  model  with  respect  to  the 
fitted  model.  Let  v{z  |  T,  N)  be  the  true  intensity  and  let  z/(  z  |  T.  N)  denote  one  of  the  fitted 
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models,  i.e., 


N 

u(z  I  T,  TV)  =  ^  AT(  z  I  h(Tj),  Sw  )  (37) 

i= 1 
N 

is(z  I  T,  N)  =  M  (  z  |  h(Tj),  £w  )  .  (38) 

1=1 

Thus,  the  true  distribution  of  z  can  be  written  as 

1  N 

Mz)  =  Sw)’  (39) 

i= 1 

and  its  approximation 


1  if  ^ 

p(z)  =  7^5^(zlh(Ti),  Ew)  '  (4°) 

P*  i=  1 

Now  the  Kullback-Leibler  information  of  p( z)  with  respect  to  p(z)  is 


Dkl(p\P)  =  /  p(z)  logp(z)  dz-  /  p(z)  logp(z)  dz, 


>T 


<T 


(41) 


which  is  a  measure  of  the  divergence  of  p( z)  relative  to  p(z).  The  aim  is  make  the 
Kullback-Leibler  information  small.  As  the  first  term  on  the  right-hand  side  of  equation  41  does 
not  depend  on  the  model,  only  the  second  term  is  relevant.  It  can  be  expressed  as 


'/({zi<  zm}\ P) 


it 


l°gp(z)  p(z)dz, 
logp(z)  dP(z), 


JT 

where  P  denotes  the  true  distribution  and  (zi,  . . . ,  zm }  is  the  observed  data, 
estimate  of  ?/({zi,  ....  zm}  \P)  is  given  by 


(42) 

(43) 

Now  the  sample 


^  IIL 

V({zi,  •  •  • ,  zm}  |P)  =  —  V  log/3(z  )  (44) 

m  ' 

3= 1 

=  —  log£  ({zh  . . . ,  zm}  |  T ,  n)  ,  (45) 

where  £(•)  indicate  the  likelihood  function  and  equation  45  follows  from  the  assumption  that  the 
observations  are  independent.  The  bias  ofrj({z1,  . . . ,  zm}  | P)  as  an  estimator  of 
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rj({ zlr  . . . ,  zm }  | P)  is  the  functional 

b(P)  —  EP  rj({zu  . ..,  zm}\P)  -77({zi,  zm}  |P)  ,  (46) 

where  EP  denotes  expectation  with  respect  to  the  true  distribution  P.  Thus,  an  information 
criterion  for  model  section  can  be  based  on  the  bias-corrected  log  likelihood  given  by 

log£  ({Zl,  . . . ,  zm}  |  T,  iv)  —  b(P),  (47) 

using  an  appropriate  estimate  of  the  bias  term.  The  intent  is  to  select  the  model  that  would 
maximize  the  above  quantity  and  thus  minimizes  the  Kullback-Leibler  information.  In  the 
literature,  the  information  criteria  formed  are  generally  expressed  in  terms  of  twice  the  negative 
value  of  the  difference  given  above,  so  they  are  of  the  form 

-21og£  ({Zl,  . . . ,  zm}  |  T,  iV)  +  2 b(P),  (48) 

The  intent  therefore  is  to  choose  a  model  to  minimize  the  criterion  equation  48. 

3.2.1  Akaike’s  Information  Criterion  (AIC) 

Akaike  (29)  showed  that  b(P)  is  asymptotically  equal  to  d,  where  d  is  equal  to  the  total  number  of 
parameters  in  the  model.  Thus  from  equation  48,  AIC  selects  the  model  that  minimizes 

—2  log  £  ({z!,  . . . ,  zmj  |  T,  iV)  +  2d.  (49) 

For  the  2-D  localization  problem  under  consideration,  the  above  criterion  can  be  rewritten  as 

-21og£({zi,  ...,  zm}  |T,  iV)  +4iV.  (50) 

Therefore,  the  present  approach  would  select  several  different  candidates  for  N  and  solve  the 
corresponding  localization  problem  using  the  EM  algorithm.  Afterwards,  the  AIC  associated 
with  each  model  is  calculated  and  the  model  corresponding  to  the  lowest  AIC  is  selected  as  the 
estimated  model.  A  summary  of  the  proposed  algorithm  is  presented  in  table  1. 
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Table  1.  Summary  of  algorithm. 

Given  data:  {zi,  . . . ,  zm} 

N 

Fit  the  intensity:  u(z  |  T,  N)  =  -A/”  ( z  I  h(Tj),  £w  ) 

2=1 

Estimate  the  parameters:  N  and  {T  | .  . . . ,  Ty  } 

•  Select  the  candidate  models'.  N\. . . . .  Nf< 

•  FOR  k  =  1  :  I< 

-  FOR i=  l  :  Nk 

*  Initialize  T ( k )  €  T 

-  END  FOR 

-  FOR  EM  iteration  n  =  0, 1, . . .  until  converges 

*  FOR  i  =  1  :  Nk 

■  Update  rj”+1)  ( k )  using  equation  35 

■  Update  Ty™+1^  (. k )  using  equation  36 

*  END  FOR 

-  END  FOR 

^  m  ^  ^ 

-  Calculate  the  maximum  likelihood:  =  exp  {— Nk}  \  [  K^|T(fc),  Nk ) 

3=1 

-  Calculate  AIC(k )  =  —  21og£fc  +  4 Nk 

•  END  FOR 

•  Select  the  index  k  corresponding  to  the  minimum  AIC: 

j  =  min  {AIC(1), . . .  ,  AIC(k), . . . ,  AIC(K)} 

k 

•  Parameter  estimates  are  N:)  and  T(j) 
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4.  Numerical  Example 


Consider  a  scenario,  where  there  are  four  targets  (. N  =  4)  located  at 


Following  the  setup  presented  in  section  2.,  the  measurement  intensity  is  given  as 

4 

u(z  I  T)  =  yjtAf(z\ h(Tj),  Ew)  +  ud(  z), 

i=  1 

[9  0 

where  the  known  measurement  covariance  matrix  is  £w  =  ,  i.e.,  the  range 

0  5  x  10”3J 

standard  deviation  amounts  to  3  m  and  the  bearing  standard  deviation  to  4°.  The  /,  ’s  are  given  as 

J1;4  =  {0.5732, 1.7962, 1.2172,  0.4135} 


and  the  clutter  intensity  is  given  as 


ud(z)  =  0.0122A f 


z  I  h 


281. 25l  \ 
0.75  )  ’ 


900  0 

0  0.4874 


The  sensor  location  is  selected  as  Sx  —  10  and  Sy  —  0.  In  order  to  evaluate  the  performance  of 
the  algorithm,  200  MC  runs  were  conducted.  For  each  run,  the  number  of  measurements  are 
determined  by  sampling  from  the  Poisson  distribution: 

(4)m 

PM\m)  «  — -  exp  {-4}  . 
ml 

The  m  measurements  are  obtained  from  the  mixture  pdf: 


4.0122 


^/?;A7(z|h(T,),  £w)  +  0.0122AT  z|h 


281.25 


0  0.4874 


with  mixing  probabilities  {0.1429,  0.4477,  0.3034,  0.1031,  0.0030}.  Here  we  use  the  optimal 
subpattem  assignment  (OSPA)  metric  to  assess  the  performance  of  the  localization 
algorithm  (30).  To  be  specific,  we  use  the  OSPA  metric  of  order  1  with  cut-off  value  100. 
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Given  in  figure  1  is  the  histogram  obtained  from  the  MC  runs  for  the  number  of  measurements 
and  estimated  number  of  targets.  Note  that  the  histogram  for  the  number  of  measurements 
clearly  resembles  that  of  a  Poisson  pdf  while  the  histogram  for  the  estimated  number  of  targets 
clearly  favors  4. 


(a)  Number  of  measurements  (b)  Estimated  number  of  targets 

Figure  1 .  Histogram  for  number  of  measurements  and  estimated  number  of  targets. 


Figure  2  contains  the  OSPA  metric  obtained  for  each  MC  runs.  Note  that  the  OSPA  metric 
around  the  100  mark  indicates  the  MC  runs  with  a  cardinality  error  of  one  while  the  OSPA  metric 
around  the  200  mark  indicates  the  MC  runs  with  a  cardinality  error  of  two.  Note  that  the 
majority  of  the  metrics  are  well  below  25  indicating  no  cardinality  error  and  accurate  localization. 


Figure  2.  OSPA  metric  for  each  MC  runs. 


Given  in  figure  3  is  the  localization  results  obtained  for  one  of  the  MC  trials.  For  the  numerical 
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results  presented  in  figure  3,  m  =  7  and  the  measurement  are 


z  = 


295.4439 

248.0845 

266.9659 

308.8356 

-7.4361 

5 

7.9683 

5 

-19.7274 

1 

17.9654 

301.3977 

251.7503 

251.5474 

-2.3297 

5 

6.9274 

5 

3.2172 

20 

15 

10 

5 

C/2 

£  o 

-5 

-10 

-15 

-20 

-100  -50  0  50  100  150  200  250  300  350 

x-axis 

Figure  3.  Single-sensor  simulation  scenario. 

For  all  MC  runs,  five  different  models  are  selected  and  their  corresponding  AIC  are  given  in 
table  2. 


O  Sensor 

♦ 

★  Target 

'  t . t  v ' ' ; . ! . ; . . 

<§>  Measurement 

Estimate 

. . : . : . . n  . : . . . 

'★ 

;  ;  ; 

i _ i _ i _  A _ i _ 

Table  2. 

Models  and  AIC. 

AIC 

N\  =  2 

A4  =  3 

Ai  =  4 

IVi  =  5 

JVi  =6 

26.7675 

24.7321 

24.3130 

27.1980 

31.8099 

Thus,  the  parameter  estimates  are  N  =  4  and 


T 


/  I” 298.4321 
\  -4.9096 


250.4694 

6.0468 


266.9659 

-19.7274 


308.8356 

17.9654 
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5.  Multi-Sensor  Problem  Formulation 


This  section  considers  a  multi-sensor  scenario  where  there  are  ‘s’  sensors  located  at  points 


/  q  q  q  N  /  S*1  S X2  S*s 

V^l)  ■  ■  ■  i  ^ s )  I  „  )  „  )  •  •  •  )  „ 

V  '-5 yi  ^V2 


The  sensor  locations  are  known  a  priori.  Following  the  problem  formulation  given  in  section  2., 
the  number  of  measurements  and  the  measurement  set  are  jointly  modeled  as  a  PPP.  Thus, 
consider  the  s  different  realizations  of  s  PPPs,  vpl1),  \p(2);  _ . . ;  \p(s): 


^(1)  =  (jn i,  {z^, 


(i)  „(i) 


Z2  > 


^ ]  =  (jnt,  {zS£),  zf,  •  •  • ,  z(g})  , 

V>(s)  =  (ms,{zSs),  z(2\  . ..,  zjj}j  . 

For  £  —  1 , ,s,  the  f-th  realization,  i.e.,  'ipU),  corresponds  to  the  number  of  measurements  and 
the  measurement  set  obtained  by  the  f-th  sensor.  The  measurement  noise  is  assumed  to  be 
zero-mean  Gaussian  with  known  variance  £w.  Now  following  the  similar  formulation  presented 
for  the  single-sensor  scenario,  the  f-th  PPP  is  modeled  via  the  intensity  function 

N 

V®  (z)  =  Y  PD  (T<) M  (z I  h  ( T* ,  Se ) ,  £w)  +  V®  (z) . 


Let  ^i:s  denotes  the  set  of  all  measurements,  i.e., 


=  Y  me,  {z™,  ..., 


(!)  7(i)  7W  7W  7W  7(s)  t  1 

1  . . z1  . ,  ■  ■  ■  ,Z1  ,  ...  ,  Zma  f  I  . 


Jmi  ?  *  *  *  ?  "1  )  *  *  *  i  ^m£  ?  *  *  *  ?  "1  ?  *  *  *  ?  "to. 


Since  the  superposition  of  two  PPPs  is  also  a  PPP,  can  be  considered  as  a  PPP  with  intensity 


s  /  N 


7  (z)  =  £  =  £  £  pD(Tl)AT(z|h(Ti,^),£w)  +  ^)(z) 


f=i  \  i=i 
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Now  p  (V,i:s)  can  be  written  as 


V  {ll>l:a)  =  exp 


IT 


s  mu 


7(z)  dz 

e= i  j=i 


M' 


(53) 


6.  Multi-Sensor  Solution 


The  multi-sensor  solution  proposed  in  this  section  directly  follows  from  the  single-sensor  solution 
presented  previously. 

6.1  EM  Method 

Given  the  number  of  targets,  N,  the  intensity  i/f)  (z  |  T)  of  the  f-th  sensor  is  the  superposition  of 
N  Gaussian  components  along  with  the  clutter  intensities: 

N 

iy{e\ z  |  T)  =  PD( Ti)  J\f(z\h  ( Tj,  Se  ) ,  sw)  +  vf  (z).  (54) 

i= 1 


6.1.1  E  -  Step 


First  consider  the  f-th  sensor,  where  there  are  rn/:  measurements.  For  j  =  1, . . . ,  me,  the  natural 
choice  of  the  “missing  data”  for  the  C-th  PPP  are  the  conditionally  independent  random  indices 
kj£\  G  (0,1,2,...,  A^},  that  identify  which  of  the  Gaussian  components,  i.e.,  target, 
generated  the  measurement  z  ‘  J .  Here  kj  =  0  indicates  that  the  measurement  is  generated  by 
clutter.  Let 

ZP  =  (rne  ,  { (zf : 1 ,  kP ),...,  (z f ,  kf ) , . . . ,  (z  ^ ) })  (55) 

denote  the  complete  data.  For  i  e  (0, 1,  2, . . . ,  N},  let 

ZP(*)  =  {(zf,*f  )  :  =  *}■  (56) 

Let  (i)  >  0  denote  the  number  of  indices  j  such  that  k(-  '  =  z,  and  let 

^P(i)  =  (nP(j),zP(i))  .  (57) 
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It  follows  from  the  definition  of  that  ipie\ i )  is  a  realization  of  the  PPP  whose  intensity  is 

w  /  ^(Tj:)AT(z|h(Tl,^),Sw),  ifi  e  {1,2, . .  .,N}; 

(*)  ~  \  (fit  s  f  ■  rt  (5g) 

{  ”cl  (Z)>  if  l  =  $. 

Now  only  considering  the  first  scenario,  %  e  {1,  2, . . . ,  N},  the  pdf  of  ( % )  can  be  written  as 


P  I  T)  =exp  (-  /  pD(Ti)J\f  (z\h(Ti,Se) ,  Ew)  dz 


T 


J]  pD( Ti)^'(zf|h(Ti,^),i:w) 


j:k\e)=i 


(59) 


and 


p  (ip®  (0)  |  T)  =  exp  (  -  vf  (z)  dz  )  JJ  vf  (z 


'3  > 


j:k^=tl 


Since  the  superposed  components  are  independent,  we  have 


(60) 


p  (zf  I  T  )=p  (4e\(b)  I  T)  p  (4e\l)  I  T)  p  (4e\2)  I  T)  .  ..p(4e\N)  |  T) 

-  exp  f-  I  i/w(z|T )dz)  JJ  pD (T k(t)) N  [zf  |h(TfeW,^),  £w  )  z 

\  JT  /  .  _  a\  3  3  .  .  ^ 


'j  > 


Let 


(,)/«).  \  J  pD(Tl„l)V(z<')|h(TM„,5(),  Ew),  if*f  e{l,2,...,Ar}; 

"  V3  1  1  pflzf),  ' 


if  =  0. 


(61) 


Then 


/  \  mi 

p  (z?  |  T)  =  exp  ^  i/(<)  (z|T)dzjn  (zf *  I  ) 


(62) 


The  log  likelihood  function  of  T  given  z^  is 
£ 


me 

(T  |  z?)  =  -  /  z/W(z  |  T)  dz  +  ]T  log  ^  (zf  |  T 
Jr  ,=i  v 


.W 


(63) 
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Now  note  that 


P 


/  r  \  171  £ 

(0^  |  T)  =  exp  J  v^\z  |  T)  clz  J  z/^  (zp  |  T 


(64) 


Thus  the  conditional  pdf  of  the  missing  data,  [kf  \  . . . ,  klP'j ,  can  be  written  as 
,  .  p(z^\T)  “<  /i(f)  (zf  |T  to) 

v  k{e)  I  tT  =  -A  '  =  TT 

P'  1  mJ0  ’  /  p(0M  |T)  11  ^(zWjt) 

Furthermore,  assuming  that  the  sensor  measurements  conditioned  on  the  target  locations  are 
independent  of  each  other  yields 

P  (z?\-  •  -  ,  ZA  I  T)  =  P  (ZA  I  T)  . .  .p  (z<s)  |  T) 


(65) 


n 

i=i 


/  /»  \ 

exp  j  vu)(z  |T)  dzj  (zf}  I  Tfcw 


(66) 

(67) 


and 


p  (0(1), . . . ,  0(s)  |  T)  =  p  (0«  |  T)  . .  .p  (0(s)  |  T) 


n 

e= i 


exp  I"  i/^(z  |  T)  dz  ^  fj  i/W  (zf  |  T 


(68) 

(69) 


Now  the  conditional  pdf  of  the  missing  data  for  all  sensors,  (  kp\  . . . ,  kml,  •  •  • ,  k[s\  . 
given  T  and  0(  1 00)  can  be  written  as 

P  (^l\  ■  ■  -  Ami  I  V'1'-  ■  ■  ■  ,  T)  = 


k ^ 

j  h ms  I  i 


p  (  z^ , . . . .  zis)  |  T 


p(0(O,...,0(*)|T) 

*  (zf  |  T  w) 

TT  TT _ Ai _ h_L 

tiM  ^}(zf|T) 


The  log  likelihood  function  of  T  given  z ^ , . . . ,  z^  is 


■C(T|zW,..,zW)  =  X; 


£=1  L 


/m£ 

v^\z  \T)  dz  +  ^  log  (zp  |  T 
i=i 


.(<) 


(70) 

(71) 


(72) 
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s 

Given  j(z  |  T)  =  1J<'  f  }  (z  |  T),  the  log  likelihood  function  can  be  rewritten  as 

t=  1 


«  s  rri£ 

£  (T  |  zf1, . . . ,  z^)  =-  7(z|T)  dz  +  ^T^T  log  /rW  (zf  |  T 

Jr  e=i  3=1 


.m 


Let  n  —  0,1,...  denote  the  EM  iteration  index,  and  let  the  current  feasible  value  of  T  be 
^(n)  _  . . . ,  V  Now  the  EM  auxiliary  function  is  the  conditional  expectation 


(s),T(n)‘ 


Q  (T  |  T<”>)  =  E  ?> . z,„k(1) . ^  w  [£(T|z«. 


N  N 

E  •••  E  £(T 

fc,(1)=0 


r(l) 


yK*J  y(s 

ZC  >  •  •  •  >  ZC 


’inn 


(a)  (  (a)  I  T(n) 

t1  1  Z/3  I  1 


7s)  —  a 


=i  pi.  ^ a)  (4a)  i  tw) 


Substituting  the  log  likelihood  yields 


Thus 


Q  (T  |  TW)  =  -  I  7(z  |  T)  dz+ 


N 


N 


E  -  E 

fci1)=0  *4&=0 


s  nri£ 

EE 

.1=1  3=1 


(£)  I  rp 

I  Tk^ 


/,(«)  ( 7(“)  |  rp(n) 
s  ma  P  l  |  1  (a 

TT  IT _ V _ 

J  M  M  i/w  (zy}  i  tw) 


(73) 


(74) 


(75) 


Note 


s  N 


N 


Q  (T  |  TW)  =  ]T  ]T  E  1oS  /*W  ( I  T, 


<=1  J=1  *f=0 

,,(«)  f  z(«)  I  T(n)  ^ 

^  ^4  1 


(0 


S  ji  (Zy}  |  tw)  -/r 


-  J  7(z|T)  dz 


N 


N 


E  -  E 

4X)=0  *&J=0 


S  (z(y  |  T  (o') 

TT  TT- _ v  kJ-i 

^  ^  ^  ^  (ri  |  t j 


1=1  7  =  1  ^ 


=  1 


(76) 


(77) 
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Thus  we  have 


s  mi  N 

Q( T|T<">)=]T£  £ 

e=i  j= i  fc<«=0 


N 

E 

lX1)  ...  fc(0  1.(0  ...  2.(s) 

K'l  5  —  ^ms—  v 

M) 


log  f!{e)  (zf  I  T 


i.W 


~(ri  I  rr(n) 

*  '  V 


x - 7—7- - J~k - /  7(z  I  T)  dz 

z/W  ( z  (f]  I  T(nA  -Ar 


(78) 


Note  that 


N 

E 

t.(0  ...  1.(0  1.(0  ...  J u(s)  — (7) 

K'l  5  1^mS~  V’ 


log  /r(£)  (zf  I  T 


1.(0 


(ZW|TW 


(0  (  zw  I  T(n) 


,(0 


^(zf  |  T^> 


fcW 


i/W  (  zf  |  TW 


(79) 


Thus,  Q  (T  |  T^)  can  be  written  as 


N  s  mi 

Q(  T|T<">)=EEE 

i=0  1=1  j=l 


(t)  I  rp 

z}  |  T, 


/iW  fzflTf 


z/M  (  zf  |  T» 


'T 


7(z  |  T)  dz  (80) 


Now  substituting  equation  61,  Q  (T  |  T^)  can  be  rewritten  as 


Q  (T  |  T(n)) 


P  s  me 

/  7(z|T)dz  +  EE  log 

•'T  B—  1  o—l 


=1  3=1 


N  s  mi 

EEE  log  pD(T07V(zf  \h(Ti,Se),  S, 


z=i  ?=i  j=i 


pD(Tf)V(z‘f||h(TS"),S(), 

z/W  ( z^}  I  TW 


(81) 


Recall  that 


7(z  |  T)  =  X)  pw(z  |  T)  =  £ 

*=i  r=i 


N 


E  PD(T.)  V  (z|h  ( T„  S, ) ,  Ew)  +  P<?(z) 


,  2=1 
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Thus,  Q  (T  |  T(n))  can  be  written  as 


Q  (T  |  T(n))  =  Y 


e= 1 


fTUd( z)  dz  +  Y  loS  vdi*T) 


{t),  {t)\ 

"cl  (zi  ) 


z/(ri  (  z{P  I  T» 


iV  s 

EE 

1=1  1=1 


Trig 

Y  los  pD (T>)  ^  ( zf  i  h(T-  #)>  sw ) 

i=i 


pD(T<"))JV'(zf  IhfTf.S,), 
z/(ri  f z(P  I  TW 


-y^^(TJ)Ar(z|h(T?:,E,),Ew)  dz 


Thus 


N 


Q  (T  |  TW)  =  YQ*  (T*  I  T(n))  +  Qci  (T(n)) 


i=  1 


(82) 


(83) 


where 

S 

Qi  (T,  I  TW)  = 


1=1 


Wit 

Y  log  PD(Ti)  A f  (  zf  I  h(Tj,  Se),  Ew  ) 

3= 1 


X 


pD(T‘">)JV-(zf  |h(T‘“>,S,), 
I/W  ( z^}  I  T(n) 


-^pD(Ti)3V'(z|h(T,,^),Ew)  dz 


(84) 


and 


«d  (T<">)  =  ^ 


1=1 


"d\z)dz  +  V  log  z^j(zH 


wit 


{£)/  (i)\ 
vd  (zi  ) 


'r  °  cl  y  3  '  i,W  ( z{P  |  T(") 


(85) 


This  completes  the  E-step. 


6.1.2  M  -  Step 

The  M-step  maximizes  Q  (T  |  T(nl)  over  all  feasible  T,  i.e., 

T(n+1)  =  argmax  Q  (T  |  T(n))  .  (86) 
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Assuming  there  is  no  functional  relation  between  T,  and  T;;  for  i  ^  j,  the  required  M-step 
maximum  is  found  by  maximizing  the  expressions  Qt  (T,:  |  T^f  separately.  Let 

fn+1)  =  arg  max  Q%  (T,:  |  T(n))  ,  1  <  i  <  N.  (87) 

Ti 


Before  we  further  proceed,  note  that  J\f  (  zf  |  h  ( T,  ,  S( ) ,  Ew  )  can  be  written  as 


•A/fzf  |h(Ti,5<),Sw)  = 


^/det(27rS^ 


:  exp  |  ~  (z  f  -  h  (  Tj,  Se 


xSwx  (zf -h(T. i,Se 


Thus 


log  AT  (  z  f  |  h  ( Th  S£  ) ,  Ew  )  =  -  log  ^det(2vrEw) 

-  i  (zf  -  h  ( T.,  s,  ))T  e;1  (zf  -  h  ( T„  s, ))  . 

Define  the  weight  cuf  (zf  |  T(n)?  Ew  j  as 


p^(Tf))7V(zf  \h(T<r\se), 

u(t)  (zf  |  TWj 


(88) 


The  weight  cuf  (zf  |  T(nT  Ewj  is  the  probability  that  the  point  zf  is  generated  by  the  z-th 
target  given  the  current  estimates  T99.  Now  equation  87  can  be  rewritten  as 


rp(TlH-l) 

i 


arg  max 

T  i 


E 


£=1  L j=l 


mu 


log  AT(zf  |h(Ti,^),  Ev 


UJ- 


W  (J*) 


T(n) ,  Ev 


+ 


log  pD( Ti)  d^  (zf  |  TW,  Ew)  -  jV( TO  M  (z|h  ( Tj,  St ) ,  Sw)  dz 


(89) 


Unlike  the  single-sensor  case,  here  we  do  not  consider  two  different  scenarios  based  on  pD{ Tj). 
Here  pD(Tj)  =  /,  are  assumed  to  be  scalar  quantities  that  need  to  be  estimated  along  with  T. 

Here  we  assume  that  pD(Ti)  are  scalar  quantities,  indicated  as  Jj,  that  need  to  be  estimated  along 
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with  Tj.  Thus  the  optimization  problem  in  equation  89  can  be  rewritten  as 


1 1  —  arg  max  y. 

Ti  t=i  L  j= i 

me 


m,( 

E  log  V  (  z f  I  h(T„  St),  E»  )  wf  (*f  I  T'">,  Ew)  + 


E  log  /iwf  (zf|T<">,Ew)  -  /  /W(z|h(T(,S,),Sw)  dz 
3  =  1 


mg 

E  log  V  (  zf  I  h(T„  5(),  £w  )  wf  (zf  I  T<”>,  Ew)  + 


/f+1)  =  arg  max 

4  *=i  L  j=i 

me  „ 

log  U  wf  (zf  |  TW,  Ew)  -  y  A/'(z|h(  Tj,  St ) ,  Ew)  dz 


3= 1 


From  equation  90,  T,-n+1)  satisfies  the  necessary  condition: 


E 

t= i 


3=1 


cuf  (zf  I  TW,  Ew)  Ewx  [zf  —  h  ( Tj,  St  )J  VTih  ( Tj,  Se )  - 
f  h  A/"(z|h  ( Tj,  S'* ) ,  Ew)  E"1  [z  -  h  ( Tj,  Se )]  VTth  ( Tj,  Se )  dz 


>T 


The  above  condition  can  be  rewritten  as 


i=  i 


=  0 


(  ^ 

mg 

X>f  (zf  |T<">,EW) 

zf  -  h  ( Tj,  St ) 

l 

-  3=1 

it 


Ii N  (z|h  ( Tj,  St ) ,  Sw)  [z  -  h  (  Tj,  St )]  dz 


VTih  ( Tj,  Se )  }  =0 


(90) 


(91) 


(92) 


(93) 


Based  on  the  assumption  \TM  (z|h  ( Tj,  St ) ,  Ew)  z  dz  «  h  (  Tj,  Se )  and 

frJ\f  ( z  |  h(Tj),  Ew  )  dz  &  1,  the  EM  updates  T -n+ 1 1  is  calculated  as  the  solution  to  the  equation 


me 


EE 

3=1 


U- 


(0  1.(0 


T(n)  v 


(ri 

Z  ;  = 


me 


1=  1 


EE 

3=1 


cu> 


(0 


T(n),Ew)  h(  T 


-,(n+l) 


St 


(94) 


i=i 


Due  to  the  nonlinear  nature  of  the  above  equation,  there  exist  no  closed  form  solution  to  above 
equation.  There  exist  several  ways  to  solve  the  above  equations  such  as  Newton’s  root  finding 
algorithm.  On  the  other  hand,  one  can  also  solve  equation  94  directly  using  a  search  method  or  a 
gradient  descent  method. 
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On  the  other  hand,  I-n+ ^ 


satisfies  the  necessary  condition: 


S 


E 


(95) 


Thus 


s  m,£ 

EE-f,(zfiT<",.sw) 

j(n+ i)  _  1=1  i=l _ 

V  [  J\f  (z|h  ( Tj,  ) ,  Ew)  dz 

<=i  •/r 


(96) 


6.2  AIC-Based  Model  Selection 


The  previous  subsection  provides  an  EM  approach  to  estimate  the  target  locations  given  the 
number  of  targets.  Presented  in  this  approach  is  an  information  theoretic  approach  to  select  the 
appropriate  model,  i.e.,  the  number  of  targets.  Similar  to  the  single-sensor  scenario,  here  we  use 
the  AIC  to  select  the  appropriate  number  of  targets. 

As  seen  in  the  E-step,  the  probability  p  ■, . . . ,  \  T)  can  be  written  as 

p(v>(i),...,V’i*iiT)=n 

i= i 

Since  JT  ]  (z  T)  c/z  =  N,  where  Ar  is  a  candidate  for  the  number  of  targets,  the  likelihood 
£  ^  T,  N  |  . . . ,  -060  j  can  be  written  as 


£(t,V|V><1|,....V’<‘))  =  n 


e= i 


exp 


-N 


mu 

n 

3= 1 


(98) 


Now  following  the  information  given  in  subsection  3.2.1,  the  AIC  can  be  calculated  as 

AIC  =  -21og£  (f.iV|^(1),...,^(s))  +  6iV.  (99) 

Therefore,  the  proposed  approach  would  select  several  different  candidates  for  N  and  solve  the 
corresponding  localization  problem  using  the  EM  algorithm.  Afterwards,  the  AIC  associated 
with  each  model  is  calculated  and  the  model  corresponding  to  the  lowest  AIC  is  selected  as  the 
estimated  model. 
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7.  Multi-Sensor  Example 


Consider  a  scenario,  where  there  are  four  targets  (N  =  4)  located  at 


298  250  267  310 


-5  ’  7  ’  -17  ’  18  r 


and  three  sensors  located  at 


10  \S, 


0.5  '  S, 


-15  &  SX3 


2.5  5, 


Following  the  setup  presented  in  section  5.,  the  measurement  intensity  corresponding  to  the  £— th 


sensor  rs  given  as 


i'(,)  (z  I  T)  =  £  pD(T,)  AT (z|h  ( T„  S, ) ,  Ew)  +  i/<?(z). 


For  simplifying  the  problem  here  we  assume,  (z)  is  the  same  for  all  £  e  {1,2,3}  and 

n  [9  0  1 

p  (T i)  =  I*.  The  known  measurement  covariance  matrix  is  £w  =  ,  i.e.,  the  range 


"  v  ‘  w  [0  5  x  10-3J  ’ 

standard  deviation  amounts  to  3  m  and  the  bearing  standard  deviation  amounts  to  4°.  The  /*’ s  are 


given  as 


J1:4  =  {0.5525, 1.8399, 1.2119,  0.3957} 


and  the  clutter  intensity  is  given  as 


uci( z)  =  0.0122A/’  z  |  h 


281.25  \  900  0 

0.75  )  ’  0  0.4874 


Thus  the  intensity  corresponding  to  the  PPP  consists  of  all  measurements  among  the  three  sensors 


is  given  as 


7(z)  =  EE  Ii  J\f  (z\h  ( Tj,  Sg ) ,  Ew)  +  vjp(z) 


1=1  \  i=i 


In  order  to  evaluate  the  performance  of  the  algorithm,  200  MC  runs  were  conducted.  For  each 
run,  the  number  of  measurements  for  each  sensor  is  determined  by  sampling  from  the  Poisson 
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distribution 


(4)m 

pM(rn)  «  exp  {-4}  . 

For  the  Ath  sensor,  the  rnf  measurements  are  obtained  from  the  mixture  pdf 

900  0  IV 

0  0.4874  J 

with  mixing  probabilities  {0.1377,  0.4586,  0.3021,  0.0986,  0.0030}.  After  obtaining  the 
measurements  for  all  three  sensors,  the  proposed  EM  algorithm  is  executed  for  each  model  using 
the  appropriate  initial  condition  selected  from  the  following  set: 


Here  we  consider  the  following  five  models: 

IVi  =  2,  N2  =  3,  N3  —  4,  N4  —  5,  and  N5  =  6. 


320 

250 

295 

265 

275 

300 

20 

5 

0 

1 

0 

5 

-10 

5 

5 

5 

15 

1 


4.0122 


M  (z|h  ( Ti?  St ) ,  Sw)  +  0.0122A/’ 


i=  1 


z|h 


281.25 

0.75 


For  each  model  the  AIC  is  then  calculated  according  to  the  information  given  in  subsection  6.2 
and  the  model  corresponding  to  the  minimum  AIC  is  selected  as  the  correct  model. 

Here  we  use  the  OSPA  metric  to  assess  the  performance  and  accuracy  of  the  proposed  localization 
algorithm  (30).  To  be  specific,  we  use  the  OSPA  metric  of  order  1  with  cut-off  value  100. 

Given  in  figure  4  is  the  histogram  obtained  from  the  MC  runs  for  the  number  of  measurements  for 
all  three  sensors.  Note  that  the  histogram  for  the  number  of  measurements  resembles  that  of  a 
Poisson  pdf  with  a  rate  of  four. 
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(c)  Sensor  3 


Figure  4.  Histogram  for  number  of  measurements. 


Figure  5a  contains  the  histogram  obtained  from  the  MC  runs  for  the  number  of  estimated  targets. 
The  histogram  for  the  estimated  number  of  targets  clearly  indicates  that  there  are  four  targets. 
Figure  5b  contains  the  OSPA  metric  obtained  for  each  MC  runs.  Note  that  the  OSPA  metric 
around  the  100  mark  indicates  the  MC  runs  with  a  cardinality  error  of  one  while  the  OSPA  metric 
around  the  200  mark  indicates  the  MC  runs  with  a  cardinality  error  of  two.  Note  that  most  of  the 
metrics  are  well  below  50  indicating  no  cardinality  error  and  accurate  localization. 


29 


(a)  Estimated  number  of  targets 


(b)  OSPA  metric  for  each  MC  runs 


Figure  5.  Estimated  number  of  targets  and  OSPA  metric  for  each  MC  runs. 


Given  in  figure  6  is  the  localization  results  obtained  for  one  of  the  MC  trials.  For  the  numerical 
results  presented  in  figure  6,  mi  =  6,  m2  =  2,  m3  =  5  and  the  measurement  are 


286.5969  6.2309 
240.8335  0.0081 
259.9712  6.2273 
302.0675  0.0577 
256.8199  6.2335 
288.3582  6.2747 


331.4559  0.0351 
310.6632  6.2480 


313.1043 

0.0071 

267.7569 

6.2744 

283.0506 

6.1852 

321.8741 

0.0963 

262.4796 

0.0491 
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(a)  Sensor  1 


(b)  Sensor  2 


★ 

Target 

o 

Sensor 

♦ 

Measure 

☆ 

Estimate 

10 


-10 


❖ 


'  it 


♦ 


-100  -50  0  50  100  150  200  250  300  350 

x-axis 


(c)  Sensor  3 


Figure  6.  Individual  sensor  measurements  and  EM  solutions. 


For  all  MC  runs,  five  different  models  are  selected  and  the  correct  model  is  selected  as  the  one 
with  lowest  AIC  score.  Thus  the  parameter  estimates  are  N  =  4  and 


T 


297.1037 

-7.3038 


311.0966 

20.7981 


267.6813 

-17.5251 


250.3221 

5.1730 


The  fused  solution  along  with  the  individual  measurements  are  given  in  figure  7,  where  the  fused 
solutions  are  displayed  as  green  stars.  As  shown  in  figure  6,  sensor  one  yields  six  measurements 
(indicated  as  blue  diamonds),  sensor  two  yields  two  measurements  (indicated  as  cyan  squares), 
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and  sensor  three  yields  four  measurements  (indicated  as  magenta  triangles). 
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Figure  7.  Multi-sensor  fused  solution. 


8.  Experimental  Results 


This  section  presents  the  implementation  of  the  proposed  algorithm  for  multi-shooter  localization 
using  a  network  of  acoustic  GDSs.  The  individual  GDSs  are  composed  of  a  passive  array  of 
microphones  that  is  able  to  localize  a  gunfire  event  by  measuring  the  direction  of  arrival  for  both 
the  acoustic  wave  generated  by  the  muzzle  blast  and  the  shockwave  generated  by  the  supersonic 
bullet  (2-5).  After  detecting  a  gunfire,  the  individual  sensors  report  their  solution,  usually  in  the 
form  of  range  and  bearing  to  the  shooter  locations  relative  to  the  sensor,  along  with  their 
orientation  and  GPS  positions  to  a  central  node  over  a  communication  network.  At  the  central 
node,  the  individual  sensor  solutions  are  fused  along  with  the  GPS  positions  to  yield  a  highly 
accurate,  geo-rectified  solution.  More  details  on  shooter  localization  using  a  network  of  acoustic 
GDSs  can  be  found  in  references  (31-37). 

Experiments  were  conducted  for  the  quad  symmetric  sensor  formation  given  in  figure  8  using  a 
sensor  network  composed  of  nine  sensors.  The  sensor  pattern  spreads  over  25  m  front  to  back. 

In  figure  8,  the  shooter  position  is  marked  by  a  red  human  figure,  and  the  shot  line  is  marked  by  a 
translucent  red  line.  The  GPS  locations  of  the  three  shooter  positions  are  given  in  table  3. 
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Figure  8.  Quad  symmetric  sensor  formation. 


Table  3.  Shooter  Locations. 


Shooter  Position 

GPS  -  East  (m) 

GPS  -  North  (m) 

Shooter  Position  1 

283309 

4709539 

Shooter  Position  2 

283270 

4709567 

Shooter  Position  3 

283337 

4709632 

The  sensor  locations  and  headings  correspond  to  the  quad  symmetric  formation  are  given  in 
table  4. 


Table  4.  Sensor  locations  and  heading  for  quad  symmetric 
formation. 


Sensor 

GPS  -  East  (m) 

GPS  -  North  (m) 

Heading  (deg) 

SW1 

283130 

4709427 

40 

SW2 

283129 

4709434 

39 

SW3 

283165 

4709401 

31 

UGS1 

283133 

4709431 

39 

UGS2 

283169 

4709398 

30 

UGS3 

283168 

4709405 

31 

VM1 

283127 

4709431 

40 

VM2 

283172 

4709402 

30 

VM3 

283177 

4709395 

29 

Here,  30  shots  were  fired  for  each  shooter  position.  For  each  run,  the  number  of  measurements 
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reported  by  each  sensor  is  modeled  as  the  Poisson  distribution 


pM(m)  « 


(3)m 


m\ 


exp  {-3}. 


Following  the  formulation  presented  in  section  5.,  the  measurements  obtained  for  the  7-th 
sensor  are  modeled  as  i.i.d.  samples  from  the  mixture  pdf 
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with  mixing  probabilities  {0.4378,  0.4853,  0.0728,  0.0041}.  The  covariance  associated  with 
individual  measurements  are  obtained  from  the  confidence  weights  provided  by  the  sensors  (32). 
After  obtaining  the  measurements  for  all  nine  sensors,  the  proposed  EM  algorithm  is  executed  for 
each  model  using  the  appropriate  initial  condition  selected.  Here  we  consider  the  following  four 
models: 

Ni  =  2,N2  =  3,  N3  =  4,  and  NA  =  5. 

For  each  model  the  AIC  is  then  calculated  according  to  the  information  given  in  subsection  6.2 
and  the  model  corresponding  to  the  minimum  AIC  is  selected  as  the  correct  model. 


Similar  to  the  previous  results,  here  also  we  used  the  OSPA  metric  of  order  1  with  cut-off  value 
100  to  assess  the  performance  and  accuracy  of  the  proposed  localization  algorithm  (30). 

Given  in  figure  9  is  the  histogram  obtained  from  the  MC  runs  for  the  number  of  measurements  for 
all  nine  sensors.  Note  that  the  histogram  for  the  number  of  measurements  resembles  that  of  a 
Poisson  pdf  with  a  rate  of  three. 
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(a)  Sensor  1 


(b)  Sensor  2 


(c)  Sensor  3 


(d)  Sensor  4 


(e)  Sensor  5 


(f)  Sensor  6 


(g)  Sensor  7  (h)  Sensor  8  (i)  Sensor  9 


Figure  9.  Histogram  for  number  of  measurements  for  each  experimental  runs. 
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Figure  10  contains  the  histogram  obtained  from  the  30  experimental  runs  for  the  number  of 
estimated  targets.  The  histogram  for  the  estimated  number  of  targets  clearly  indicates  that  there 
are  three  targets. 
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Figure  10.  Estimated  number  of  targets. 


Figure  1 1  contains  the  OSPA  metric  obtained  for  each  experimental  runs.  Note  that  most  of  the 
metrics  are  well  below  20  indicating  no  cardinality  error  and  very  accurate  localization. 
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Figure  1 1 .  OSPA  metric  for  experimental  runs. 


Given  in  figure  12  is  the  localization  results  obtained  for  one  of  the  experimental  trials.  For  the 
numerical  results  presented  in  figure  6,  m\  =  2,  m2  =  6,  m3  =  1,  m4  =  4,  m5  =  4,  m6  =  2, 
m7  =  4,  m8  =  3,  and  m9  =  3.  In  figure  12,  the  measurements  obtained  for  each  sensors  are 
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denoted  as  blue  diamonds,  the  true  shooter  locations  are  denoted  as  red  stars,  the  sensor  location 
is  denoted  as  yellow  circles,  and  the  estimated  shooter  locations  are  represented  as  green  stars. 
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Figure  12.  Experimental  results  for  individual  sensor  measurements  and  EM  solutions. 


The  fused  solution  along  with  the  individual  measurements  are  given  in  figure  13,  where  the  fused 
solutions  are  displayed  as  green  stars.  As  shown  in  figure  13,  the  fused  solution  is  more  accurate 
compared  to  the  individual  measurements. 
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Figure  13.  Experimental  result  for  multi-sensor  fusion. 


9.  Conclusion 


This  report  presents  the  finite  point  process  approach  to  the  multi-target  localization  problem  for 
the  single-sensor  as  well  as  multi-sensor  scenario.  Here  it  is  assumed  that  target  identification  is 
not  possible,  and  therefore,  no  association  between  the  measurements  and  targets  are  available. 
Furthermore,  the  number  of  targets  in  the  surveillance  region  is  unknown,  and  due  to  the  limited 
range  of  the  sensors,  missed  detections  can  occur  and  the  presence  of  clutter  can  induce  false 
alarms.  Here  we  propose  an  EM  algorithm  to  estimate  the  target  locations  while  the  information 
criterion,  AIC,  is  used  to  estimate  the  number  of  targets.  The  preliminary  implementation  of  the 
proposed  algorithm  on  synthetic  data  produced  accurate  results.  Implementation  of  the  proposed 
algorithm  on  experimental  data  obtained  for  the  multi-shooter  localization  problem  further 
confirms  the  numerical  results.  Here  we  use  the  optimal  subpattem  assignment  metric  of  order  1 
with  a  cut-off  value  of  100  to  assess  the  performance  of  the  localization  algorithm.  As  the  results 
given  in  sections  7.  and  8.  indicate,  the  finite  point  process  approach  is  able  to  accurately  estimate 
the  number  of  targets  and  their  locations  in  the  presence  of  clutter  and  missed  detection.  Future 
work  will  include  decreasing  the  sensitivity  of  the  proposed  algorithm  to  the  initial  guess  with  the 
use  of  multiple-thread  search  and  deterministic  annealing  EM  algorithm  (38).  The  current 
scheme  can  also  benefit  from  a  systematic  approach  for  selecting  potential  target  models. 
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